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Abstract 

The criticality of the (2 + 1) -dimensional 5* = 1 transverse-field Ising 
model is investigated with the numerical diagonalization method. The scaling 
behavior is improved by tuning the coupling-constant parameters; the S = 1 
spin model allows us to incorporate a variety of interactions such as the single- 
ion anisotropy and the biquadratic interactions. Simulating the clusters with 
= 6, 8, . . . , 20 spins, we estimate the critical indices as 1/z/ = 1.586(6) and 
T] = 0.034(4). 

Keywords: 

75.10.Jm 75.40.Mg 05.50. +q , 05.70.Jk 

1. Introduction 

The finite-size-scaling method allows us to estimate the critical indices 
from the simulation data. In practice, however, there appear sub-leading 
singularities, the so-called corrections to scaling. Corrections to scaling are 
suppressed gradually by enlarging the system sizes. On the one hand, there 
was reported a scheme to control corrections to scaling directly by tuning 
the coupling-constant parameters [H, ^ [s], 0, Is], S, 0, S M 10 1- Particularly, 



in the lattice-field theory, the lattice artifact brings about a severe obstacle 
to extracting information through the continuum limit. There have been 
reported a number of attempts to suppress the lattice artifact by tuning ex- 
tended interaction parameters; see Ref. (sj for a review. This issue is closely 
related to the analysis of criticality, because the critical-point phenomena 
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are also governed by long-wave-length fluctuations. In practice, it is even 
desirable to develop a scheme where no ad hoc adjustment of interac- 
tion parameters is required. This approach may be particularly of use for 
the numerical diagonalization method, with which the tractable system size 
is restricted severely. (On the contrary, as for the Monte Carlo method, it 
might be more rewarding to enlarge the system size rather than incorporate 
extra interactions.) 



Motivated by the above idea, in Ref. [Uj, we investigated the three- 
dimensional (classical) Ising model with the extended (next-nearest-neighbor 
and plaquette-four-spin) interactions by means of the diagonalization (transfer- 
matrix) method for the clusters with < 13 spins. Thereby, we estimated 
the critical indices as z/ = 0.6245(28) and i/h = 2.4709(73). In the present 
paper, based on this preliminary survey, we shall make the following improve- 
ments. First, we consider the quantized version of the Ising model, namely, 
the (2 + l)-dimensional transverse-field Ising model. (Note that the (2 + 1)- 
dimensional ground-state phase transition belongs to the same universality 
class as that of the three-dimensional classical counterpart.) The quantum 
Hamiltonian matrix is sparse (the matrix has few non-zero elements), as com- 
pared to the transfer matrix of the classical model. Taking this advantage, 
we treat the system sizes up to < 20. Second, we extend the magnitude 
of constituent spin to S = 1 from S = 1/2. Owing to this extension, we are 
able to incorporate a variety of interactions such as the single-ion anisotropy 
D and the biquadratic interactions, J4 and J4. To be specific, we consider 
the Hamiltonian 

n = -jj^s^s^-j' sts^-.hY,{sts^^f-j',Y,{sts^^Y+DY^{stf-TY,si 

(1) 

Here, the quantum S = 1 operators {Sj} are placed at each square-lattice 
point i. The summations, ^^^^ possible nearest- 

neighbor and next-nearest-neighbor pairs, respectively. The parameters J 
and J' are the corresponding coupling constants. The parameter F denotes 
the transverse magnetic field. Last, these coupling-constant parameters are 
adjusted to 

{J, J' , Ji, J'^, D) = (0.41191697085,0.16125069616,-0.11764020018, 

-0.05267926601, -0.39781956122). (2) 

This finely adjusted parameter set renders a significant improvement to the 
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scaling behavior. We determine the parameter set in Sec. [21 based on a 
real-space decimation (Fig. [T]). The F-driven phase transition is our concern 
(Sec. ED. 

In fairness, it has to be mentioned that the Ising model with 5 = 1 and 



D = \n2 was investigated [l2j by means of the Monte Carlo method. At 



D = ln2, the cluster-update algorithm works efficiently. We stress that in 
our simulation, the anisotropy D is tuned so as to suppress the finite-size 
errors; technically, the value of D is not subjected to any restriction in the 
numerical-diagonalization scheme. 

The rest of this paper is organized as follows. In Sec. [2], setting up a real- 
space-decimation scheme, we determine a (non-trivial) fixed point ([2]). In 
Sec. [31 around the fixed point, we perform a large-scale computer simulation; 
the simulation scheme is explicated in the Appendix. Taking the advantage 
of suppressed scaling errors, we analyze the criticality of the d = 3 Ising 
universality class with the diagonalization method. In Sec. [H we preset the 
summary and discussions, referring to related studies. 



2. Search for a scale-invariant point: Suppression of finite-size 
corrections 

As mentioned in the Introduction, the coupling-constant parameters, 
( J, J', J4, J4, D), are set to Eq. ([2|). In this section, we explain how this 
parameter set is determined. Here, we stress that the analysis (real-space dec- 
imation) of this section is not intended to determine the critical exponents. 
The critical exponents are determined in the subsequent scaling analysis in 
Sec. m 

To begin with, we outline the underlying idea of the finely-tuned param- 
eters, Eq. ([2|). The parameter set ([2|) is a scale-invariant (fixed) point of 
a real-space decimation. In Fig. [H we present a schematic drawing of the 
real-space decimation for a couple of clusters labeled by 5* and L. The cluster 
sizes are 2x2 and 4x4, respectively. As indicated, through a decimation of 
the • spins, the L cluster reduces to a coarse-grained one identical to the S 
cluster. Our aim is to search for a scale-invariant point with respect to the 
real-space decimation. 

Before going into the fixed-point analysis, we set up the simulation scheme 
for the clusters S and L. We cast the Hamiltonian into the following plaquette- 
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based expression 



i 1 



IX 



(3) 



[ijkl] i i 



with the plaquette interaction 



■□ 

ijkl 



( Q-^ _L _L _L Q'^\ 1' ( _l 

J',{{StStf + {S^^Slf). 



(4) 



(The denominator 2 compensates the duphcated sum.) Hence, the Hamilto- 
nian for the S cluster is given by 



Here, the parameter h controls the boundary-interaction strength, and here- 
after, we set h = 0.75; we consider the validity of this choice afterward. The 
boundary-interaction parameter b interpolates smoothly the open {b = 0) 
and periodic (6 = 1) boundary conditions. The point is that for such a small 
cluster with the linear dimension L = 2, the bulk and boundary interactions 
are indistinguishable. Such a redundancy is absorbed into the tunable param- 
eter b] in other worlds, the parameter b is freely tunable without violating the 
translation invariance. We make use of this redundancy to obtain a reliable 
fixed-point location. On the other hand the L cluster does not have such a 
redundancy, and the Hamiltonian is given by Eq. ([T]) with L = 4 unambigu- 
ously. We diagonalize these Hamiltonian matrices numerically. Note that 
we employ the conventional diagonalization method, rather than Novotny's 
diagonalization method, which is utilized in Sec. |3l see the Appendix as well. 

With use of the simulation technique developed above, we search for the 
fixed point of the real-space decimation. We survey the parameter space, 
regarding F as a unit of energy; namely, we set F = 1 throughout this 
section. Thereby, we impose the following conditions 



4 4 



'Hs = {i + &)?/° 34 + D J2ist? -rJ2s!. 



(5) 



i=l i=l 



{SJSJ)l 



(6) 
(7) 
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{{sts^r)s = {is^Ar)L (8) 

{{SlSir)s = {{Sl~Slf)L (9) 
{{Slf)s = {{Slf)L. (10) 

as a scale-invariance criterion. The arrangement of the spin variables is 
shown in Fig. [H the symbol (. . .)s,l denotes the ground-state average for the 
respective clusters. The above equations set the scale-invariance condition as 



to the correlation functions [13|. In order to solve the non-linear equations 
O-dinD, we utilized the Newton method, and arrived at the fixed point, 
Eq. ([2]); the last digits may be uncertain because of the round off errors. 
In the next section, we present the simulation result around the fixed point 
([2]). As claimed in Ref. |14l], the decimation scheme based on the correlation 
function(s) is problematic. In our approach, the quantitative analyses of 
critical exponents are made subsequently by the finite-size scaling. 

Last, we argue the validity of the fixed point ([2]), and the boundary 
condition h = 0.75. In Sec. 13. we estimate Fc = 1.0007(17) via the finite- 
size-scaling analysis. Clearly, this result is consistent with F = 1 postulated 
in this section (decimation analysis). As a matter of fact, the above real-space 
decimation contains several sources of errors, such as the restricted system 



size and a fundamental difficulty as to the determination of fixed point [15 
Here, we make use of the redundant parameter b, aiming to control these 
errors in a rather ad hoc manner. 



3. Numerical results 

In this section, we present the simulation result for the (2 + l)-dimensional 
Ising model ([T]). The coupling-constant parameters are set to the scaling- 
invariant point (|2]), at which the scaling behavior improves. (We stress that 
the scale- invariant analysis of Sec. [2] is simply a preliminary one, and it is no 
longer used.) We employ Novotny's numerical diagonalization method; the 
technical details are addressed in the Appendix. Owing to this method, we 
treat a variety of system sizes = 6, 8, . . . , 20 (A^ is the number of spins); 
note that conventionally the system size is restricted within A^ = 4, 9, 16, . . .. 
The linear dimension L of the cluster is given by 

L = y/N, (11) 

because the A^ spins constitute a rectangular cluster. 
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Conventionally, the d = 3 criticality has been studied by means of the 
Monte Carlo method (rather than the diagonalization method). The error of 
the Monte Carlo simulation comes from a purely statistical origin. On the 
one hand, the diagonalization result is free from such a statistical error, and 
the data analysis is not straightforward. In the present paper, we perform 
two independent analyses for each critical exponent, aiming to appreciate 
possible systematic errors. 

3.1. Critical point r c 

In this section, we provide evidence of the F-driven phase transition with 
the other coupling constants set to Eq. 

In Fig. [21 we plot the scaled energy gap LAE for various F, and = 
6,8,..., 20. According to the finite-size-scaling theory, the scaled energy gap 
LAE should be invariant at the critical point. Indeed, we observe an onset 
of the F-driven phase transition around F ~ 1. 

In Fig. [31 we plot the approximate transition point Fc(Li, L2) for [2/ + 
L2)f with 6 < A^i < A/'s < 20 and Li,2 = a/A^2; the vahdity of the 
extrapolation scheme (abscissa scale) is considered at the end of this section. 
The approximate transition point Fc(-Li,-L2) denotes a scale-invariant point 
with respect to a pair of system sizes (Li, L2). Namely, according to the phe- 



nomenological renormalization 16|], the approximate transition point satisfies 
the equation 

LiAE(Li)|r=r4L„L.) = L2AE(L2)|r=r.(L„L.). (12) 

The least-squares fit to the data of Fig. [3] yields Fc = 1.0007(17) in the 
thermodynamic limit L — )■ cxd. As to this estimate, F^ f[T2|) . possible system- 
atic errors are not appreciated, and the amount of error margin is spurious. 
Because the estimate Fc is not used in the subsequent analyses of critical 
exponents, we do not go into further details; rather, Tc{Li,Ls) is used. 

We mention a few remarks. First, we consider the abscissa scale, 
utilized in Fig. [3l Naively, the scaling theory predicts that the dominant 



corrections to Fc should scale like l/L'^+^Z'' 17| with u = 0.821(5) and 1/u = 
1.5868(3) [13]. On the one hand, in our simulation, such dominant corrections 
should be suppressed by adjusting the coupling constants to Eq. ([2l). The 
convergence to the thermodynamic limit may be accelerated. (For extremely 
large system sizes, the singularity l/L'^'^^^'^ may emerge.) Hence, in Fig. [3l 
we set the abscissa scale to Second, we argue a consistency between 
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the finite-size scaling and the real-space decimation (Sec. [2]). In Sec. [21 we 
made a fixed-point analysis, regarding F as a unit of energy; namely, we set 
r = 1. This proposition is quite consistent with the above simulation result, 
Tc = 1.0007(17), justifying the fixed-point analysis in Sec. El In other words, 
around the fixed point, corrections to scaling may cancel out satisfactorily. 
Third, the shaky character of Fig. [31 is an artifact due to the cluster size; 
as the cluster size deviates from the quadratic number (A^ = 9, 16,25), the 
shaky character becomes amplified. This feature is intrinsic to Novotny's 
diagonalization method. (A bump of Fig. [31 should be attributed to this 
character.) 

3.2. Critical exponent rj 

In Sec. 13. 1[ we observed an onset of the F-driven phase transition. In this 
section, we calculate the critical exponent i], based on the finite-size scaling. 
In Fig. [H we plot the approximate critical exponent 

rj{Li,L2) = -21n[m(Li)/m(L2)]|r=r.(Li,L,)/ln(Li/L2) - 1, (13) 

for [2/(Li + La)]^ with 6 < A^i < A^'a < 20 (Li,2 = ^A^); afterward, we 
consider the abscissa scale, 1/L^. Here, the magnetization m is given by 
m = {{g\M'^\g)y^'^ /N with the ground state \g) and the total magnetization 
M = S^. The least-squares fit to these data yields i] = 0.03088(26) in 
the thermodynamic limit L — )• oo. 

In Fig. [5l we plot the approximate critical exponent 

r]{L,, L2) = - ln[x(Li)/x(L2)]|r=r.(L„L.)/ ln(Li, L^) + 2, (14) 

for [2/(Li + ^2)]^ with 6 < A^i < A''2 < 20. Here, the susceptibility x is given 
by the resolvent form 

with the ground state energy Eg. The resolvent form is readily calculated 



with use of the continued-fraction method [19[. The least-squares fit to the 
data in Fig. [51 yields 77 = 0.03676(19) in the thermodynamic limit L — )• 00. 
Considering the difference between the results of Figs. [H and \5\ as an error 
indicator, we estimate the critical exponent as 

rj = 0.034(4). (16) 
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The index immediately yields the critical exponent 



P/v = 0.517(2), (17) 

through the scaling relation. 

We address a number of remarks. First, we consider the error margin of 
the critical exponent. The result t] = 0.03676(19) (Fig. |5]) is in remarkable 
agreement with that of recent Monte Carlo, i] = 0.0368(2) [l2|. In order 
to appreciate the error margin properly, one observation would not suffice, 
because there should exist an error (systematic deviation) that cannot be 
captured by a statistical (least-squares fit) analysis. Therefore, we take into 
account another independent simulation result rj = 0.03088(26) (Fig. H]). (In 
general, the reliability of A^-independent observations is given by a/y/N — 1 
with the standard deviation a of each observation.) That is, the deviation 
of the results Figs. H] and [5] yields an indicator of the error margin. Hence, 
we arrive at the above conclusion with the error margin, which covers these 
independent observations. Second, we argue the abscissa scale 1/L^ utilized 



in Figs, mandm In Ref. [18| (see, in particular. Fig. 3 and Sec. IIB), it was 
claimed that the abscissa scale 1/L^ yields a consistent result; namely, results 
via two independent analyses coincide. Hence, we set the abscissa scale to 
1 /L^ in Figs. H] and El Naively, the scaling theory predicts that the leading 
corrections to the critical exponent should scale like with u = 0.821(5) 
(lij . On the other hand, the convergence should be accelerated more than 
the naively expected one, because such dominant corrections are truncated 
by adjusting the coupling-constant parameters to Eq. (|2]). In this sense, an 
accelerated convergence 1/L^ is a reasonable one, at least, within the range 
of system sizes tractable with the diagonalization method. 

3.3. Critical exponent v 

In this section, we calculate the critical exponent v. 
In Fig. ini we plot the approximate critical exponent 

ln(^e(Li)/(9re(L2))|r=r,(Li,L2) , 3 

21n(L,/L,) + 2' ^'^^ 

for [2/(Li -|- -^2)]^ with % < Nx < < 20. Here, the internal energy, e, is 
given by the off-diagonal part of the Hamiltonian (namely, e = —{g\^ ^j^i ^i\9) 1-^)^ 
which turns out to exhibit little size dependence. The least-squares fit to 
these data yields 1/v = 1.58765(25) in the thermodynamic limit L 00. 
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In order to obtain the exponent 1/z/, we have to take a F derivative with 
respect to a certain quantity. In the preliminary stage, we surveyed a vari- 
ety of quantities such as the internal energy set tentatively to e = [diagonal 
part of Eq. ([I])]/iV, excitation gap, and susceptibility. As a consequence, we 
arrived at the conclusion that the above-mentioned clue f lTSj) exhibits little 



size dependence. Similar tendency is observed [2l| for the case of = 1/2, 
namely, the conventional transverse- field Ising model. In the case of 5 = 1/2, 
via the same scheme, we obtained a tentative result 1/z/ 1.6 even with- 
out making any extrapolation. This fact indicates that the F derivative 
dre = Co + AqL?'/"^^ has a negligible constant term Co, and the singular 
term Aq is dominating. In the analysis of Fig. El we made use of this char- 
acter. 

In Fig. [71 we plot the approximate critical exponent 

1/z/ - 13 /u = ln(arm(Li)/9rm(L2)|r=r.(Li,L,))/ln(Li/L2), (19) 

for [2/{Li + L2)f with Q < Ni < N2 < 20. The least-squares fit to these data 
yields 1/z/ — /3/z/ = 1.0676(55) in the thermodynamic limit L — )■ 00. Com- 
bining this result with Eq. f lTTl) . we arrive at 1/z/ = 1.5846(59). Because it is 
calculated indirectly, it might have uncontrolled deviation. Nevertheless, the 
consistency to the above result (Fig. Ej) seems to be satisfactory. Counting 
these independent results, we estimate 

1/z/ = 1.586(6), (20) 

with the error margin coming from the latter one. 

3.4. Simulation at (J, J', J4, J^, D) = (0.6, 0, 0, 0, 0) 

In the above, we simulated the transverse- field Ising model, Eq. ([T|), with 
the finely tuned coupling constants, Eq. As a comparison, in this section, 
we provide a result for a model without the extended interactions; namely, 
we set (J, J, J, J, D) = (0.6, 0, 0, 0, 0) tentatively. 

In Fig. [HI we plot the scaled energy gap LAE for various F and = 
6, 8, ... , 20. We observe an onset of the F-driven phase transition around 
F ~ 1. However, the location of critical point appears to be unclear, as 
compared to that of Fig. [21 This result demonstrates clearly that the finely- 
tuned coupling constants lead to an elimination of finite-size corrections. 



9 



4. Summary and discussions 



The (2 + l)-dimensional S = 1 transverse-field Ising model ([T]) was inves- 
tigated numerically. Unlike the conventional S = 1/2 spin model, the S = 1 
model allows us to incorporate a variety of interactions, (J, J', J4, J4, D). By 
adjusting these interactions to a scale-invariant point ([2]), we attain substan- 
tial elimination of scaling errors; see Figs. [2] and [HI As a result, we estimate 
the critical indices as 77 = 0.034(4) and 1/p = 1.586(6). Through the scaling 
relations, we arrive at 



We overview related studies. As mentioned in the Introduction, a set 



of indices [0.1265(84), 0.3304(48), 1.213(11)] was obtained in our prelimi- 



nary survey [ll| as for the Ising model with the next-nearest-neighbor and 



plaquette- four-spin interactions. The extension of interactions, Eq. ([T]), 
leads to an improvement as to the indices (12T]) . Moreover, in the present 
study, we manage two independent analyses of each index. Hence, it is 
expected that the error margins of Eq. (12T]) would be appreciated prop- 



erly. A numerical diagonalization result [0.1144(24), 0.3281(13), 1.2319(65)] 



|20| was evaluated for a large 6x6 cluster (without the extended inter- 



actions) for 5* = 1/2; a diagonalization of a 7 x 7 cluster [N = 49) far 
exceeds a limitation of available computer resources. The present scheme 
(see the Appendix) admits us to enlarge linearly (for example, A^ = 
6,8,10,...), and an extension of A^ might not be so computationally de- 
manding. The series-expansion method yields the following indices via the 



high-temperature expansion [0.1096(5), 0.32653(10), 1.2373(2)] [2^ as well 
as the field-theoretical analyses, [0.109(4), 0.3258(14), 1.2396(13)] H and 
[0.1091(24), 0.3257(5), 1.2403(8)] 0. Recent Monte Carlo results are [0.1109 
(15), 0.3262(4), 1.2366(15)] Q and [0.10940(36), 0.326695(88), 1.23721(27)] 



[I2I ]. The latter result, which has considerably small error margins, was esti- 
mated by surveying a variety of lattice structures. Actually, the character of 
scaling corrections depends on the respective lattice structures. Hence, it is 
sensible to count various lattices in order to reduce and certify the extrapola- 
tion errors; in the present paper, we carry out two independent analyses for 
each critical exponent. At present, such a posteriori consideration cannot be 
omitted. In this sense, further consideration of scaling corrections would be 
desirable in order to estimate the critical indices in a fully controlled manner. 



(a,/3,7) 



[0.1084(72), 0.3260(18), 1.2396(53)]. 



(21) 






10 



L 



Figure 1; A schematic drawing of the real-space renormahzation group (decimation) for 
the (2 + l)-dimensional transverse-field Ising model ([T|) is presented. Through decimating 
out the spin variables indicated by the symbol • within the L cluster, we obtain a coarse- 
grained lattice identical to the S cluster. Imposing the scale-invariance conditions, Eqs. 
(|6))-([T0|. we arrive at the fixed point ([2]). 



3.5 I 1 1 1 1 1 r 




0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 

r 

Figure 2: Scaled energy gap LA_E is plotted for various D and N — 6, 8, . . . , 20 
[L = \/N); note that we survey the F-driven phase transition with the other interac- 
tions (J, J', J4, Ji,D) adjusted to a fixed point ([2]). We observe a clear indication of the 
F-driven transition around F ^ 1. Clearly, the scaling behavior is improved as compared 
to that of the conventional transverse-field Ising model without the extended interactions 
(Fig. ©. 
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Figure 3: The approximate critical point Fc ([T^ is plotted for [2/{Li + i2)]^ with 6 < 
Ni<N2< 20 (ii^2 = ^/Nia). The least-squares fit to these data yields F^ = 1.0007(17) 
in the thermodynamic limit L — >■ oo. 
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Figure 4: The approximate critical exponent ri{Li,L2) is plotted for [2/(Li + ^2)]^ 
with 6 < TVi < iV2 < 20 (^1,2 = ^^^^1,2)- The least-squares fit to these data yields 
T] = 0.03088(26) in the thermodynamic limit L — > 00. 
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Figure 5: The approximate critical exponent ri{Li,L2) (fT4l) is plotted for [2/{Li + ^2)]^ 
with 6 < iVi < A^2 < 20 (ii,2 = y^Vi^). The least-squares fit to these data yields 
r] = 0.03676(19) in the thermodynamic limit L — ^ 00. 
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Figure 6: The approximate critical exponent Xjv ([18]) is plotted for [2/(Li + ^2)]^ with 
6 < iVi < A^2 < 20 (ii,2 = ^A^i,2)- The least-squares fit to these data yields Xjv = 
1.58765(25) in the thermodynamic limit L — >■ 00. 



13 



an. 
> 




1.05 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 

n2 



[2/(Li+L2)]^ 



Figure 7: The approximate critical exponent 1/v — P/v (IT^ is plotted for [2/(Li + ^2)]^ 
with 6 < iVi < iV2 < 20 (-Zji,2 = The least-squares fit to these data yields 

1/v ~ P/v — f .0676(55) in the thermodynamic limit L — )■ 00. 




Figure 8: Tentatively, we turned off the extended interactions, {J, J' , J4, J^^, D) = 
(0.6, 0, 0, 0, 0), and calculated the scaled energy gap LAE for various T and TV = 6, 8, . . . , 20 
{L — Vn). We notice that the data are scattered as compared to those of Fig. [2l 
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Appendix A. Numerical method: Novotny's diagonalization scheme 



In Sec. 131 we simulated the two-dimensional transverse- field Ising model 

In this Appendix, we 



([T]) with Novotny's diagonalization method [26|, 127 , 
explain the simulation algorithm. 

To begin with, we outline the simulation algorithm. The spins constitute 
a one-dimensional [d = 1) alignment rather than a. d = 2 cluster. The 
dimensionality is lifted to = 2 effectively by introducing the long-range 
interactions with the (average) distance \/N. Because of this geometrical 
character, one is able to construct a cluster with an arbitrary number of spins 
N = 6,8,..., 20; note that conventionally, the system sizes are restricted 
within A^ = 4,9,16,.... 

Originally, Novotny's method was developed for the classical Ising model. 
Meanwhile, it was extended (formulated) to adopt the transverse magnetic 



field [28|; see the Appendix of Ref. [29|] as well. In the present paper, we 
resort to the formalism addressed in Ref. [2^ . In the following, we summarize 
the (almost trivial) modifications in order to extend the spin magnitude 
to S* = 1 from S = 1/2. First, the notation for the spin operator, af'^'^ 
(Pauli matrices), has to be replaced with S'f Correspondingly, the Hilbert- 
space basis [Eq. (Al) of Ref. 
with i = 1,2, ... ,N and Si = 
[summand of Eq. (A5) in Ref. 



29j] is represented by the expression K^j}) 
—1,0,1. Second, the plaquette interaction 
has to be replaced with an extended one 
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{{s^}\Tm) 



^ T 
- J:((5,T,+i)' + iSj+iTjY) 



(A.l) 



Technically, the above modifications suffice for the simulation of the 
Hamiltonian ([T]). In the following, we provide a supplementary remark con- 
cerning the translational operator P; see Eq. (A4) of Ref. 29|. We calculate 
the matrix element of the operator P'" with v = VN, based on the expansion 



(A.2) 



2TTkv 



{{S^}\Plm) = J2J2^{S.}\'^,)exp{^^){^,m), 
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with the intermediate states and an integer index k specifying a Brillouin 
zone (reciprocal space); note that the reciprocal space {k} is suitable for 
representing P. The point is that each Brillouin zone is no longer equivalent 
because of the incommensurate oscillating factor exp{2TTkvi/N) in Eq. f lA.2p . 
The range of Brillouin zone of the present simulation is —N/2 < k < N/2. 
More specifically, the intermediate sum, ^fc, in Eq. f lA.2p is expanded into 

flfc = -jO,-N/2 + (^-N/2+l + 0.-NI2+2 + • • • + aAr/2-1 + ^(i'Nl2-, (A. 3) 

k 

for a summand Ofc. The factor 1/2 compensates the duplicated sum at the 
zone boundaries. 
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